## Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan

library(R2jags)

J <- 4 ## number of issue areas asked about
K <- 5 ## number of endorsement groups (including control)
N <- 1500 ## number of survey respondents
n.province <- 4
n.per.province <- rep(N/n.province,n.province)
province <- rep(1:n.province,N)[1:N]
b.province <- c(-1,-0.25,0.25,1)
theta.province <- matrix(c(rep(0,n.province), c(0.3,0.3,0.3,0.3), c(0,0.3,0.3,0.8), c(-0.8,-0.8,0,0.8), c(-0.8,-0.3,-0.3,0)), n.province, K)
sigma <- matrix(c(rep(0,n.province), c(0.2,0.2,0.2,0.2), c(0.3,0.3,0.3,0.3), c(0.2,0.4,0.2,0.4), c(0.2,0.3,0.2,0.3)), n.province, K)


## Ideal Points
x <- rep(NA, N)
for (i in 1:N){ x[i] <- rnorm(1, mean=b.province[province[i]], sd=1) }
x <- x / sd(x)

beta <- runif(J, 0.5, 2) ## discrimination parameter
alpha <- matrix(0.1,J,4)
for (j in 1:J){  ## variation in thresholds, but distance ensured
	while(alpha[j,2]-alpha[j,1] < 0.5 | alpha[j,3]-alpha[j,2] < 0.5 | alpha[j,4]-alpha[j,3] < 0.5){
		tmp <- rnorm(4, mean=0, sd=1)
		alpha[j,] <- tmp[order(tmp)]
	}}

lambda.province <- array(NA, c(n.province, J, K))
for (i in 1:n.province){
	for (j in 1:J){
		for (k in 1:K){
			lambda.province[i,j,k] <- rnorm(1, mean=theta.province[i,k], sd=sigma[i,k])
		}}}

omega <- matrix(c(rep(0,J), runif(J*(K-1), 0.3, 1)), J, K)

s <- array(NA, c(N,J,K)) ##Individual level impact of group endorsement k
for (i in 1:N){
	for (j in 1:J){
		for (k in 1:K){
			s[i,j,k] <- rnorm(1, mean=lambda.province[province[i],j,k], sd=omega[j,k])
		}}}


##Ordered Responses (5 choices)
Y <- array(NA, c(N,J,K))
for (i in 1:N){
	treat <- ifelse(sample(c(0,1),1)==1,TRUE,FALSE) ## Half of Jake's respondents receive the control, remaining half split among treatments
	for (j in 1:J){
		for (k in 1:K){
			p1 <- pnorm(alpha[j,1] - beta[j]*(x[i] + s[i,j,k]))
			p2 <- pnorm(alpha[j,2] - beta[j]*(x[i] + s[i,j,k])) - (p1)
			p3 <- pnorm(alpha[j,3] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2)
			p4 <- pnorm(alpha[j,4] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2 + p3)
			p5 <- 1 - (p1 + p2 + p3 + p4)
			
			p <- c(p1,p2,p3,p4,p5)
			Y[i,j,k] <- sample(c(1:5), size=1, prob=p) ## Full array of potential outcomes
		}}
	if (treat==FALSE) { Y[i,,2:K] <- NA}
	if (treat==TRUE) {
		Y[i,,1] <- NA 
		for (j in 1:J){ 
			answered <- sample(2:K,1); Y[i,j,-answered] <- NA  
		}
	}}

Y2 <- matrix(NA,N,J)
endorser <- matrix(NA, N,J)
for (i in 1:N){
	for (j in 1:J){
		endorser[i,j] <- (1:5)[!is.na(Y[i,j,])]
		Y2[i,j] <- Y[i,j,endorser[i,j]]
	}}
Y <- Y2



lambda.province.real <- lambda.province

#######################################
##  Prep and send data through JAGS  ##
#######################################
data <- list ("N", "J", "K", "Y", "endorser", "province", "n.province", "lambda.province.real") 
inits <- function() {list(alpha0=matrix(seq(-2,2,length.out=(4*J)),4,J), beta=runif(J, 0.2, 2), 
			x=rnorm(N), s=matrix(rep(0, N*J), N,J), 
			delta.province=rep(0,n.province),
			lambda.province=array(c(rep(0,n.province*J),rep(0,n.province*J*(K-1))), c(n.province, J, K)), 
			theta.province=matrix(c(rep(0,n.province),rep(0,n.province*(K-1))),n.province,K),
			omega=matrix(c(rep(0,J),runif(J*(K-1))),J,K) )} 
params <- c("lambda.province.error") 
model <- "ModelSimulationsIssueProvince.txt" 
fit <- jags(data, inits, params, n.iter=20000, model.file=model)
fit <- autojags(fit, n.iter=20000, n.thin=20, Rhat=1.02, n.update=5)
save.image(file = Outfile)

